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Abstract 

Discretization  and  linearization  of  the  steady-state  Navier-Stokes  equations  gives  rise  to  a 
nonsymmetric  indefinite  linear  system  of  equations.  In  this  paper,  we  introduce  preconditioning 
techniques  for  such  systems  with  the  property  that  the  eigenvalues  of  the  preconditioned  ma¬ 
trices  are  bounded  independently  of  the  mesh  size  used  in  the  discretization.  We  confirm  and 
supplement  these  analytic  results  with  a  series  of  numerical  experiments  indicating  that  Krylov 
subspace  iterative  methods  for  nonsymmetric  systems  display  rates  of  convergence  that  are  in¬ 
dependent  of  the  mesh  parameter.  In  addition,  we  show  that  preconditioning  costs  can  be  kept 
small  by  using  iterative  methods  for  some  intermediate  steps  performed  by  the  preconditioner. 
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1.  Introduction.  Consider  the  steady-state  Navier- Stokes  problem:  given  data  f,  find 
the  velocity  u  and  pressure  p  satisfying 

—  z/V2u  -| — u(div  u)  +  u  •  Vu  +  gradp  =  f 

(1.1)  2  in 

div  u  =  0 

subject  to  boundary  conditions  on  dfl]  fl  C  R2  or  fi  C  R3.  Here,  the  scalar  v  is  the  inverse 
of  the  Reynolds  number,  or  the  ratio  of  convection  to  diffusion  in  the  system.  In  the  diffusion 
dominated  case  (//  — ►  oo)  (1.1)  tends  to  a  linear  self-adjoint  system  of  equations — the  Stokes 
problem. 

There  are  two  ways  of  calculating  solutions  to  the  system  (1.1).  A  popular  approach  is  to 
compute  “true”  steady-state  solutions  of  the  time-dependent  Navier-Stokes  equations.  There 
are  many  ways  to  do  this:  one  way  is  to  make  use  of  the  “characteristics”  associated  with  the 
hyperbolic  part  of  the  Navier-Stokes  operator  via  a  Lagrange- Galerkin  approach  (for  example, 
see  [12]).  The  associated  transpose-diffusion  splitting  leads  to  absolutely  stable  temporal  dis¬ 
cretizations  so  that  large  time  steps  can  be  taken.  At  each  time  step,  a  symmetric  indefinite 
matrix  system  corresponding  to  a  time-discretized  Stokes-like  system  must  be  solved.  These 
systems  can  be  be  solved  efficiently  by  iterative  methods,  for  example,  if  a  multigrid  solver 
is  used  to  precondition  the  primary  (Laplacian)  operator.  There  are,  however,  a  number  of 
disadvantages  to  the  time-dependent  approach.  Simple  time  discretization  methods  based  on 
the  /2-projection  onto  the  discretely  divergence-free  subspace  [9]  have  an  0(h)  CFL  restriction 
on  the  time  step,  which  impinges  on  efficiency.  On  the  other  hand,  absolutely  stable  schemes 
like  the  method  of  backward  characteristics  are  known  to  be  sensitive  to  implementation  issues 
(e.g.,  the  need  to  perform  quadrature,  see  [12]).  Even  with  fixed  grids,  efficiency  is  often  limited 
by  the  costs  associated  with  interpolation. 

In  this  work,  we  consider  the  alternative  approach  of  attacking  the  system  (1.1)  directly. 
Applying  a  fixed  point  (or  Picard)  iteration,  the  system  (1.1)  reduces  to  solving  a  sequence 
of  linear  Oseen  problems  of  the  form:  given  some  (divergence-free)  velocity  field  w,  find  the 
velocity  u  and  pressure  p  satisfying 

—  z/V2u  +  w  •  Vu  +  gradp  =  f 

(1.2)  in 

div  u  =  0 

subject  to  the  same  boundary  conditions. 

For  this  methodology  to  be  effective  it  is  necessary  to  solve  the  discrete  versions  of  (1.2) 
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efficiently.  Thus,  our  general  starting  point  is  the  matrix  problem 


(1.3) 


uA  +  N  Bf  \  f  u\  _  f  f\ 

b  o  )  Vpy  “  v ^ / 


where  A  =  Af  represents  diffusion  (for  example,  —V2),  and  hence  is  a  positive  definite  matrix 
of  order  n„,  N  represents  convection  (w  -V),  and  the  np  X  nu  matrix  B  represents  the  coupling 
between  the  discrete  velocity  u  and  the  pressure  p.  Note  that  the  representation  of  the  quadratic 
convection  term  in  (1.1)  ensures  that  N  =  —  Nf,  that  is,  the  discrete  form  of  the  convection 
operator  is  skew-symmetric  [9],  p.  53.  Note  that  if  normal  velocities  are  specified  everywhere 
on  the  boundary  then  the  system  (1.3)  is  singular,  pressure  is  only  unique  up  to  a  (hydrostatic) 
constant.  We  assume  in  the  following  analysis  that  the  pressure  solution  is  uniquely  specified 
in  this  case,  e.g.,  by  insisting  that  its  mean  is  zero. 

Working  in  a  conventional  mixed  finite  element  framework,  we  will  further  assume  that  the 
underlying  velocity  and  pressure  approximations  are  (div-)stahle  (see  e.g.,  [2],  p.  57,  [9],  pp. 
lOff,  [18]),  i.e. ,  defining  a  mesh  parameter  h,  a  velocity  space  and  a  pressure  space  Ph ,  there 
exist  constants  7,  T,  independent  of  h,  such  that 


2  <  (P,BA  1Btp)  <  r2 

ip,  Qp) 


Vp  E  Ph 


Here,  Q  is  the  pressure  mass  matrix,  or  alternatively  the  Grammian  matrix  of  basis  functions 
defining  Ph .  The  lower  bound  7  is  the  so-called  inf-sup  constant.  The  relation  (1.4)  is  cru¬ 
cial  to  the  success  of  iterative  solvers  for  solving  discrete  Stokes  problems  for  it  implies  that, 
using  a  quasi-uniform  mesh,  the  Schur  complement  BA~1Bt  has  condition  number  bounded 
independently  of  h.  It  is  also  known  from  our  previous  work  [15]  that  when  v  00,  “opti¬ 
mal”  preconditioners  for  the  Laplacian  sub-blocks  give  rise  to  “optimal”  preconditioners  for  the 
Stokes  problem  in  the  sense  that  the  spectra  of  the  underlying  discrete  operators  are  contained 
in  small  clusters,  which  are  bounded  independently  of  h.  A  consequence  of  this  is  that  the 
asymptotic  rate  of  convergence  of  Krylov  subspace  methods  applied  to  discrete  Stokes  problems 
is  also  independent  of  h. 

In  this  paper  we  derive  analogous  results  in  the  general  Oseen  case.  We  introduce  two 
preconditioners  for  the  Oseen  problem  such  that,  for  any  value  0  <  v  <  00,  the  eigenvalues 
of  the  preconditioned  Oseen  operator  are  bounded  independently  of  the  mesh  size.  These 
observations  apply  to  arbitrary  discretizations  satisfying  (1.4).  In  addition,  we  show  in  a  series 
of  numerical  experiments  that  these  bounds  on  eigenvalues  are  predictive  of  the  performance  of 
Krylov  subspace  iterative  methods  for  solving  the  preconditioned  Oseen  equations.  Of  course, 
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it  is  well  known  that  when  convection  dominates  (i.e. ,  when  v  is  “small”  relative  to  h  and 
1 1 w ||),  the  standard  Galerkin  approximation  deteriorates.  Oscillations  in  the  discrete  velocity 
are  apparent  if  the  local  mesh  Reynolds  number  Reh  =  /j||w||/zz  is  greater  than  unity.  In 
such  situations,  the  addition  of  stream  wise  diffusion  to  the  discrete  system  is  known  to  give 
added  stability,  both  theoretically  and  numerically,  see  [3]  and  [11].  In  our  experiments,  we 
demonstrate  the  effectiveness  of  the  ideas  using  both  a  standard  Galerkin  discretization  on  a 
set  of  quasi-uniform  grids,  and  a  streamline-upwind  scheme  on  a  set  of  uniform  grids. 

The  remainder  of  the  paper  is  divided  into  three  sections.  Our  main  theoretical  results 
are  presented  in  Section  2,  and  results  of  numerical  experiments  confirming  and  augmenting 
the  theoretical  analysis  are  given  in  Section  3.  In  Section  4,  we  consider  more  practical  precon¬ 
ditioning  strategies  and  present  a  perturbation  analysis  and  additional  numerical  experiments 
demonstrating  their  effectiveness. 

2.  Preconditioning  strategies.  In  this  section,  we  introduce  two  preconditioning  tech¬ 
niques  for  (1.1)  and  present  an  analysis  showing  that  the  spectra  of  the  preconditioned  systems 
are  bounded  independently  of  the  discretization  mesh  size  h.  Throughout  the  section,  we  will 
be  concerned  with  the  eigenvalues  of  preconditioned  matrices;  these  matrices  can  be  viewed  as 
being  of  the  form  AM~X  where  A  is  the  original  matrix  and  Ai  is  the  preconditioner.  Equiv¬ 
alently,  we  are  concerned  with  the  solution  of  the  generalized  eigenvalue  problem  Av  =  XAiv. 
All  the  matrices  in  question  are  implicitly  parameterized  by  h.  For  simplicity,  we  state  our 
results  under  the  assumption  that  B  of  (1.3)  has  full  rank. 

The  first  idea  is  derived  from  a  method  developed  in  [13,  15,  19]  for  the  discrete  Stokes 
equations,  where  the  coefficient  matrix  has  the  form 


(2.1) 


A  Bt 
B  0 


Consider  the  preconditioner 

A  0 
0  Q 

for  (2.1).  The  eigenvalues  of  the  preconditioned  operator  are  then  given  by  the  solution  to  the 
generalized  eigenvalue  problem 


(2.2) 


A  B<\  (  u\  _  f  A  0 
B  0  )\p)-X\ 0  Q 


One  solution  is  A  =  1,  of  multiplicity  nu  —  np,  for  which  the  eigenvectors  have  the  form  J 
where  Bu  =  0,  i.e.,  u  is  “discretely  divergence  free.”  The  remaining  eigenvalues  come  from  the 
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solution  of  the  quadratic  equation  A  (A  —  1)  =  p,  where  ^  is  a  generalized  eigenvalue  of  the 
Schur  complement  associated  with  (2.1), 

(2.3)  BA~1Btp  =  pQp. 


Equivalently, 

(2.4) 


A  = 


1  A  \/l  +  4  p 

2 


Since  (1.4)  implies  that  as  h  — ►  0,  the  solutions  to  (2.3)  remain  bounded  above  and  below,  it 
follows  that  the  eigenvalues  of  (2.2)  are  also  bounded.  The  preconditioned  conjugate  residual 
method  can  then  be  used  to  solve  (2.1),  with  a  convergence  rate  independent  of  h  [13,  15]. 

A  natural  generalization  for  the  discrete  Oseen  equations  uses  the  block  preconditioner 


(2.5) 


F  0 

o  iQ 


where  F  =  vA  +  N.  As  above,  the  generalized  eigenvalues  for 

<2-61  (b  o)(p)=a 

are  either  A  =  1  or  (2.4),  where  p  is  now  a  solution  to  the  generalized  eigenvalue  problem 


F 

0 


0 

zQ 


(2.7)  Sp  =  p(^Q^p, 

with  S  =  BF~1Bt ,  the  Schur  complement  for  the  discrete  Oseen  operator.  The  following  result, 
which  generalizes  the  analysis  for  the  Stokes  operator  in  [18],  provides  a  bound. 


Theorem  1.  Fhe  eigenvalues  of  the  generalized  Schur  complement  problem  (2.7)  for  the  Oseen 
operator  are  contained  in  a  rectangular  box  in  the  right  half  plane  whose  borders  are  bounded 
independently  of  h. 

Proof.  Let  C  =  B(— — ^ — )_Bt  denote  the  symmetric  part  of  A,  and  R  =  B(— — — )_Bt  its 
skew-symmetric  part,  so  that  S  =  C  +  R.  By  Bendixson’s  Theorem  ([16],  p.  418),  any  eigenvalue 
p  of  the  problem  (2.7)  satisfies 


(n  •  (p,Cp)  /  -r,,  ,  /  (p,Cp)  |T  ,  ^  \(p,Rp)\ 

(2.8)  mm  - — j-—  <  R e(p)  <  max  - — j-—  ,  Im  /i  <  max  - — j-—  . 

p  ( p,-Qp )  p  ( p,-Qp )  p  ( p,-Qp ) 

To  construct  bounds  on  these  Rayleigh  quotients,  it  will  be  convenient  to  refer  to  S oo  = 
BA~1Bt ,  the  Schur  complement  for  the  Stokes  operator.  For  the  symmetric  part  C  in  (2.8),  we 
use  the  relation 


(p,Cp)  =  (p,Cp)  (PtSqqp) 
(PiiQp)  (p^Soop)  (p,Qp) 
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In  light  of  (1.4),  we  need  only  consider  the  first  quotient  on  the  right  in  (2.9).  Note  that 


(2.10) 


F 


-l 


F 


—t 


=  F~ 


(F  +  F1 


F 


—  t 


=  (yA  +  N )_1  (zM)  (zM  -  N) 
1  -i 
v 


-l 


=  A-1/2  (z/J  -  ifV2)  'A”1/2, 


where  N  =  A  1!2N  A  1//2.  Consequently, 

(p,Cp)  =  (p,  BA-1/2  (z/J  -  liV2)-1  A-1!2 Flip)  =  (y,(l-±N2)-iv) 
(p^Soop)  (p,  ^BA~1Btp)  (v,v) 


where  v  =  A~x!2  Btp.  But  N  is  skew-symmetric,  so  that  the  eigenvalues  of  —  N2  are  real  and 
nonnegative.  Moreover,  since  N  and  A  are  first-order  and  second-order  operators,  respectively, 
the  eigenvalues  of  N  are  uniformly  bounded  in  modulus  by  a  constant  b  that  is  independent 
of  h  [5].  Therefore,  the  spectrum  of  I  —  yiN2  is  contained  in  the  interval  [l,  1  +  <52/z/2],  or, 
equivalently, 

<  ip,  Cp) 

b2  +  v2  ~  (p, ISoop)  ~ 

Combining  this  with  (1.4)  and  (2.9)  gives 

72^2  <  (p,  Cp)  2 

b2  +p2  ~  (p,  iQp)  ~ 

For  the  skew-symmetric  part  R  in  (2.8),  the  analogue  of  (2.9)  is 

(P,  Rp)  =  (P,  Rp)  (. P,SqqP ) 

(PiiQp)  (p,  is<x>p)  (p,Qp) 

and  as  in  (2.10),  we  have 

- - — —  =  -A~1/2  (z/J  +  IV)-1  N  {yl  -  N)~x  A~1/2. 


Therefore 


(2.11) 


(p,  Rp)  v(v,Nv) 

(P, \Soop)  (v,  (z/2/  —  N2)v) 


where  v  =  (vl  —  IV)-1  A-1/2 BT p.  The  skew-symmetric  matrix  N  admits  a  decomposition  of 
the  form  N  =  iUAUH  where  A  is  a  real  diagonal  matrix  and  U  is  unitary.  Consequently, 
N 2  =  —  UA2IJH ,  and  the  modulus  of  the  Rayleigh  quotient  on  the  right  side  of  (2.11)  can  be 
expressed  in  the  form 

v  |(w,  Aw)\ 

(w,  ( v2I  +  A 2)w) 
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This  is  bounded  by 


max  — - —  =  max  — - — 

-S<X<S  V2  +  A2  0<  A<5  V2  +  A2 


It  follows  from  elementary  calculus  that  this  maximum  is  1/2,  obtained  when  A  =  is,  giving 


\(p,Rp)\  K 
ip,\Qp)  ~  2 


The  following  result  follows  immediately  from  Theorem  1  and  (2.4). 

Corollary  1.  The  eigenvalues  of  the  discrete  Oseen  operator  preconditioned  by  (2.5)  consist  of 
A  =  1  of  multiplicity  nu  —  np,  together  with  four  sets  consisting  of  points  of  the  form  1  +  (a  ±  bi ) 
and  —a  ±  bi.  These  sets  can  be  enclosed  in  two  rectangular  regions  that  are  symmetric  with 
respect  to  Re  (A)  =  j,  whose  borders  are  bounded  independently  of  h. 

The  inclusion  regions  for  these  eigenvalues  consist  of  the  image  of  the  box  +  ,  T2  X  [—  ,  ^-] 

under  the  mapping  p  i— ►  Ml1)  given  by  (2.4).  It  can  be  shown  that  the  rectangular  regions  of 
this  result  are  contained  in 


X  [—1,  /]  and 


X  [~t,t] 


in  the  right  and  left  half  sides,  respectively,  of  the  complex  plane,  where 


The  fact  that  the  eigenvalues  for  the  preconditioned  system  derived  from  (2.5)  lie  on  both 
sides  of  the  imaginary  axis  is  a  potential  disadvantage  of  this  idea.  An  alternative  that  avoids 
this  problem  is  the  block  triangular  operator 


(2.12) 


F  Bt  \ 

o  -\Q)' 


For  this  choice,  the  preconditioned  eigenvalue  problem  is 


(2.13) 


F  B*\  ( u\  _  ( F  Bt 

B  0){p)-X{  0  -IQ 


Again,  one  solution  is  A  =  1,  now  of  multiplicity  nu.  If  A  1,  then  premultiplying  the  first 
block  row  of  (2.13)  by  BF -1  and  using  the  relation  Bu  =  —A (^Q)p  leads  to  the  equation  (2.7) 
for  the  other  eigenvalues.  Thus,  we  have  the  following  result. 
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Theorem  2.  The  eigenvalues  of  the  discrete  Oseen  operator  preconditioned  by  (2.12)  consist 
of  A  =  1  together  with  the  generalized  eigenvalues  of  S  in  (2.7).  Therefore,  the  eigenvalues  are 
bounded  independently  of  h. 

Remark  1.  Use  of  either  preconditioning  operator  (2.5)  or  (2.12)  entails  the  computation  of 
the  action  of  _F_1  at  each  step  of  an  iterative  procedure.  F  is  a  discrete  convection- diffusion 
operator  and  applying  _F_1  to  a  vector  using  direct  methods  will  be  expensive.  An  alternative 
is  to  replace  this  computation  with  an  approximation  obtained  by  iterative  solution  of  the 
convection- diffusion  equation.  We  will  examine  this  approach  in  Section  4. 

Remark  2.  Both  preconditioners  also  require  the  action  of  Q-1,  which  may  also  be  expensive, 
depending  on  the  choice  of  pressure  discretization.  In  this  case,  however,  it  is  known  that  Q 
can  be  replaced  by  some  approximation  Q  without  affecting  asymptotic  convergence  properties; 
only  the  constants  7  and  T  of  (1.4)  change  [20].  I11  the  experiments  discussed  in  Sections  3  and 
4,  we  replace  Q  with  a  diagonal  matrix  consisting  of  the  main  diagonal  of  Q. 

3.  Numerical  results  I:  Exact  convection-diffusion  solves.  I11  this  section,  we 
present  the  results  of  numerical  experiments  indicating  that  the  analysis  of  Section  2  is  predictive 
of  the  performance  of  iterative  methods  for  solving  (1.3).  Unless  otherwise  stated,  computations 
were  performed  using  MATLAB  4.1  011  a  SUN  Sparcstation-10. 
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Fig.  1.  Magnitude  and  direction  of  the  convecting  flow. 


Our  test  problem  is  a  “leaky”  two-dimensional  lid-driven  cavity  problem  in  a  square  domain 
(—  1  <  ,r  <  1  :  —  1  <  j/  <  1 ).  The  boundary  conditions  are  ux  =  uy  =  0  011  the  three  fixed 
walls  (x  =  —1,  y  =  —1,  x  =  1),  and  tix  =  1,  uy  =  0  on  the  moving  wall  (y  =  1).  The 
hydrostatic  pressure  is  not  explicitly  specified,  so  that  all  the  linear  equation  systems  we  solve 
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below  are  singular  with  a  one- dimensional  nullspace.  The  convective  “wind”  is  a  circular  vortex 
as  illustrated  in  Fig.  1,  and  is  given  by 

wx  =  2y(l-  x2 ) 
wy  =  —'2x(  1  -  y2). 

The  fact  that  there  is  no  dominant  flow  direction  makes  this  a  challenging  test  problem.  Note 
that  in  the  corners  and  in  the  center  of  the  flow  region  the  driving  flow  is  stagnant. 


Streamlines  :  equally  spaced  Streamlines  :  selected 


Fig.  2.  Uniform  64  X  64  grid  ://=!. 


Streamlines  :  equally  spaced  Streamlines  :  selected 


Fig.  3.  Non-uniform  64  X  64  grid  :  v  =  1/100. 


tJuless  otherwise  specified,  we  consider  three  values  of  the  viscosity  parameter  //,  namely 
1,  1/10  and  1/100.  When  //  =  1  we  have  diffusion  dominated  (essentially  Stokes)  flow,  whereas 
as  v  —  0  the  flow  becomes  dominated  by  the  “wind.”  Typical  flow  solutions  are  illustrated 
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in  Figs*  2  and  3.  Note  that  as  the  viscosity  is  decreased,  the  center  of  primary  recirculation 
moves  to  the  right  (the  Stokes  flow  solution  is  perfectly  symmetric  about  the  line  x  =  0),  and 
secondary  vortices  are  generated  in  the  two  bottom  corners. 

To  discretize  (1.2),  we  take  a  finite  element  subdivision  based  on  n  X  n  grids  of  rectangular 
elements.  Bearing  in  mind  the  nature  of  the  flow  solution  being  computed,  results  for  two 
representative  discretizations  are  presented  here:  a  conventional  Galerkin  approach  using  a 
quasi-uniform  sequence  of  grids,  and  a  streamline-upwind  method  using  uniform  grids  of  square 
elements  of  size  h  =  '2/n.  In  either  case,  the  mixed  finite  element  used  was  the  div-stable  “Taylor- 
Hood”  method  based  on  continuous  bilinear  pressure  with  a  continuous  bilinear  velocity  held 
defined  on  four  element  macro-elements  (see  e.g.,  [9],  p.  30). 

For  the  Galerkin  discretization,  the  quasi-uniform  grids  are  chosen  to  resolve  the  details  of 
the  how  in  the  four  corners  of  the  domain:  they  are  symmetric  about  x  =  0  and  y  =  0,  and 
in  each  quadrant  the  grid  lines  expand  uniformly  outwards.  The  64  X  64  grid  is  shown  in  the 
pressure  solution  plot  in  Fig.  4.  The  analytic  pressure  solution  is  singular  at  the  top  corners 
where  the  imposed  velocity  is  discontinuous. 


Fig.  4.  Pressure  solution  for  v  =  1/10. 

The  streamline-upwind  discretization  is  as  described  in  [11],  p.  185.  In  this  case,  the  block 
convection- diffusion  operator  F  is  perturbed  by  a  symmetric  positive  semi-dehnite  matrix  Aw. 
That  is,  F  =  —v +  Aw  +  N,  where  is  the  discrete  Laplacian  obtained  from  the  usual 
Galerkin  formulation.  Aw  is  the  discrete  form  of  a  stabilizing  term  a(w  •  Vu,w  •  Vv)  that 
adds  a  =  0(h)  diffusion  along  the  streamlines.  For  our  experiments  with  streamline  upwinding, 
we  took  a  =  hj 4.  Note  that  the  perturbation  does  not  affect  the  skew-symmetric  part  of  the 
convection- diffusion  operator,  so  that  the  analysis  of  Section  2  holds;  only  the  definition  of  the 
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diffusion  matrix”  A  is  changed,  from  to  — to  — +  ^Aup 


We  first  consider  the  bounds  of  Theorem  1.  Table  1  shows  the  extreme  real  parts  and  max¬ 
imum  imaginary  parts  of  the  generalized  eigenvalues  (2.7)  of  the  Scliur  complement  operator, 
for  v  =  1/10  and  1/100  with  the  streamline-upwind  discretization,  on  three  meshes.  (Eigen¬ 
values  for  the  64  X  64  grid  were  computed  on  a  SUN  630MP  using  Matlab  4.1.)  The  small 
changes  in  all  values  are  in  accordance  with  the  analysis,  although  it  appears  that  finer  meshes 
would  be  needed  to  produce  constant  values.  The  analysis  also  shows  that  the  real  parts  and 
largest  imaginary  parts  of  the  eigenvalues  are  bounded  independently  of  v\  the  bound  for  the 
smallest  real  part  is  proportional  to  v2 .  The  data  of  Table  1  are  in  agreement  with  the  upper 
bounds.  Figure  5  plots  the  smallest  real  parts  on  a  logarithmic  scale,  for  the  streamline-upwind 
discretization  on  a  64  X  64  grid  and  v  =  1/20,  1/40,  1/80,  and  1/160.  The  results  indicate  that 
the  lower  bound  is  also  tight. 


v  =  1/10 

v  =  1/100 

Grid 

Min  Re 

Max  Re 

Max  Im 

Min  Re 

Max  Re 

Max  Im 

16  X  16 

7.17E-2 

1.11 

0.46 

1.66E-2 

1.07 

0.20 

32  x  32 

8.75E-2 

1.64 

0.71 

1.33E-2 

1.11 

0.50 

64  x  64 

9.08E-2 

2.00 

0.87 

1.14E-2 

1.37 

0.74 

Table  1.  Eigenvalues  of  the  Scliur  complement,  for  streamline-upwind  discretization. 


Fig.  5.  Minimum  real  parts  of  eigenvalues  of  the  Scliur  complement,  for 
streamline-upwind  discretization  on  a  64  X  64  grid. 


To  solve  the  preconditioned  linear  systems,  we  consider  two  popular  Krylov  subspace  iter- 
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ative  methods  for  nonsymmetric  systems:  the  restarted  generalized  minimum  residual  method 
GMRES(s),  where  s  is  the  restart  length  [14],  and  a  simple  implementation  of  the  quasi¬ 
minimum  residual  method  (QMR)  [7]  based  on  coupled  two-term  recurrences  without  look¬ 
ahead.  We  use  s  =  10;  for  this  choice,  the  storage  requirements  of  restarted  GMRES  and 
(our  version  of)  QMR  are  essentially  identical.  Each  step  of  QMR  requires  twice  as  many 
matrix- vector  products  and  preconditioning  operations  as  one  step  of  GMRES(IO);  since  the 
preconditioning  costs  dominate,  QMR  has  roughly  twice  the  cost  per  iteration.  In  all  cases, 
we  use  right-oriented  preconditioning,  and  our  convergence  criterion  is  a  reduction  of  10-6  in 
the  /2-norm  of  the  residual,  starting  from  a  zero  initial  guess.  Random  initial  guesses  gave 
comparable  iteration  counts  in  all  cases. 

We  consider  the  standard  Galerkin  method  on  the  quasi-uniform  grid  sequence  first.  The 
iteration  counts  using  the  block  diagonal  preconditioner  (2.5)  are  shown  in  Table  2.  For  GM- 
RES(10),  the  residual  norm  is  computed  only  every  10  steps.  With  either  iterative  solver,  the 
counts  demonstrate  that  grid-independent  convergence  rates  are  obtainable  in  practice.  When 
v  is  0(1)  the  results  are  very  clean.  For  a  fixed  grid,  the  iteration  counts  grow  (linearly)  as  v 
tends  to  zero,  as  might  be  anticipated  from  the  analytical  bounds  of  Section  2.  When  v  =  1/100 
the  iteration  counts  slowly  increase  as  h  is  decreased,  again  suggesting  that  finer  grids  are  re¬ 
quired  to  see  the  asymptotic  behavior.  In  light  of  its  smaller  cost  per  step,  GMRES(IO)  is  more 
efficient  than  QMR  when  v  is  large;  however,  QMR  becomes  dramatically  more  efficient  when 
convection  becomes  dominant. 


GMRES(IO) 

QMR 

Grid 

16  x  16 

32  x  32 

64  x  64 

16  x  16 

32  x  32 

64  x  64 

v  =  1 

50 

50 

40 

43 

43 

41 

v  =  1/10 

90 

120 

120 

51 

68 

78 

v  =  1/100 

>  500 

>  500 

>  500 

143 

246 

375 

Table  2.  Iteration  counts  for  Galerkin  discretization  with  block  diagonal  preconditioner. 

The  same  observations  are  appropriate  in  the  case  of  the  block  triangular  preconditioner 

(2.12) ,  see  Table  3.  An  interesting  feature  here  is  that  for  both  Krylov  subspace  methods,  the 
number  of  iterations  is  roughly  halved  when  (2.12)  is  used  in  place  of  (2.5).  The  cost  per  step 
of  the  block  triangular  preconditioner  is  only  slightly  higher  than  that  of  the  block  diagonal 
preconditioner;  only  an  extra  multiplication  by  Bf  is  needed.  Thus,  the  triangular  method 

(2.12)  is  more  effective. 
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GMRES(IO) 

QMR 

Grid 

16  x  16 

32  x  32 

64  x  64 

16  x  16 

32  x  32 

64  x  64 

v  =  1 

20 

20 

20 

22 

22 

22 

v  =  1/10 

30 

40 

40 

28 

36 

39 

v  =  1/100 

320 

450 

>  500 

73 

126 

189 

Table  3.  Iteration  counts  for  Galerkin  discretization  with  block  triangular  preconditioner. 

Using  the  streamline  upwind  discretization  on  a  uniform  grid  sequence  gave  the  iteration 
counts  in  Tables  4  and  5.  The  results  are  qualitatively  similar  to  the  Galerkin  results  of  Tables 
2  and  3.  In  contrast,  for  a  “poor  discretization,”  i.e. ,  standard  Galerkin  on  the  uniform  grid 
sequence,  the  iteration  counts  tended  to  significantly  increase  if  the  local  mesh  Reynolds  number 
was  not  kept  in  check. 


GMRES(IO) 

QMR 

Grid 

16  x  16 

32  x  32 

64  x  64 

16  x  16 

32  x  32 

64  x  64 

v  =  1 

70 

60 

50 

49 

51 

47 

v  =  1/10 

100 

120 

120 

78 

91 

80 

v  =  1/100 

400 

>  500 

>  500 

154 

249 

382 

Table  4.  Iteration  counts  for  streamline-upwind  discretization  with  block 

diagonal  preconditioner. 


GMRES(IO) 

QMR 

Grid 

16  x  16 

32  x  32 

64  x  64 

16  x  16 

32  x  32 

64  x  64 

v  =  1 

30 

30 

30 

25 

27 

25 

v  =  1/10 

40 

50 

50 

36 

44 

42 

v  =  1/100 

180 

320 

470 

76 

131 

190 

Table  5.  Iteration  counts  for  streamline-upwind  discretization  with  block 

triangular  preconditioner. 

Remark  3.  In  addition  to  the  implementation  of  QMR  with  a  coupled  two-term  recurrence 
(QMR2)  discussed  above,  we  tested  a  version  without  look-ahead  based  on  a  three-term  recur¬ 
rence  (QMR3)  [6],  and  the  definitive  (Fortran)  implementation  of  two-term  QMR  with  look¬ 
ahead  (QMR2)  from  the  QMRPAK  directory  in  Netlib.  For  these  preconditioners,  the  perfor¬ 
mances  of  the  three  variants  were  virtually  identical.  However,  with  the  inexact  preconditioners 
of  the  next  section,  we  found  QMR2  to  be  much  more  robust  than  QMR3. 

4.  Numerical  results  II:  Inexact  convection-diffusion  solves.  The  dominant  costs 
of  the  preconditioners  of  Sections  2  and  3  come  from  applying  the  action  of  _F_1,  and  for 
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QMR,  -F_t,  to  some  vector  v  at  each  step  of  the  iteration.  In  this  section,  we  show  that  this 
operation  can  be  replaced  by  an  inexpensive  one  derived  from  an  approximation  to  T1-1 ,  with 
little  degradation  of  performance  of  the  Krylov  subspace  methods.  The  idea  is  to  replace  the 
preconditioning  operators  (2.5)  and  (2.12)  with 


(4.1) 


F  0 
0  lQ 


and 


F  Bt  \ 

o  -IQJ 


respectively,  where  F  ~  F.  Our  choice  of  F  will  be  implicitly  determined  by  the  use  of  iterative 
methods  to  compute  approximate  solutions  to  the  systems  Fw  =  v  and  Ftw  =  v,  although  the 
methodology  is  not  restricted  to  this  choice.  We  will  refer  to  the  preconditioners  that  use  the 
exact  action  of  T1-1  as  the  exact  versions,  and  those  based  on  approximations  as  in  (4.1)  as 
inexact  versions. 

An  analysis  of  the  effects  of  the  inexact  preconditioners  is  derived  from  matrix  perturbation 
theory.  Let  Q v  =  IQ.  The  preconditioned  matrix  for  the  exact  block  diagonal  preconditioner 
(2.5)  is 


An 


F  Bf\  f  F  0  \  1 

b  o  y  v°  QJ 


I 

BF-1  0  J  ’ 


and  that  derived  from  the  inexact  version  is 

Ad=(b  o)  (o  <3°„)  =ao  +  £d, 

where,  with  E  =  F  —  F, 

F  ef -1  o\ 

tD  yBF-'EF-1  0 )’ 

Similarly,  the  preconditioned  matrices  for  the  exact  and  inexact  block  tridiagonal  precondition¬ 
ers  (2.12)  satisfy  At  =  At  +  A,  where 


(I  0  \  _  (  EF-1  EF^BtQ-1  \ 

y BF-1  BF-1  B^-1  J  ’  T  yBF^EF-1  BF-1  EF-1  BtQ~1  )  ' 


We  have  the  following  bounds  on  the  eigenvalues  of  the  preconditioned  systems  using  inexact 
preconditioners. 

Theorem  3.  If  An  =  VdA^V))1  is  diagonalizable,  then  for  any  eigenvalue  fi  £  a(An), 


min  |A  —  /i|  <  \\EF  1  HooKooCA)  max(l, \\BF  1  ||oo)- 
A ) 

If  At  =  VrAyVy1  is  diagonalizable,  then  for  any  eigenvalue  ji  £  <t(At), 


min  |A-/i|<|| EF  1  ||ooKoo(Vt)  (1  +  H-AQ  1  (A)  max(l,  || BF  ^A). 
Ae  cr(AT) 
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Proof.  The  result  is  an  immediate  consequence  of  the  Bauer-Fike  Theorem  [8],  p.  342,  which 
states  that  for  diagonalizable  A  =  VAV-1,  any  g  £  a(A-\-£)  satisfies  min^e<7(^)  | A  —  /Li |  < 
k(V)  ||£||,  where  ||  •  ||  is  any  /p-norm.  □ 

Thus,  if  F  is  a  good  enough  approximation  to  F,  i.e. ,  if  enough  inner  iterations  are  used, 
then  ||ifi7’_1||  will  be  small  and  the  eigenvalues  of  Ad  and  At  will  be  close  to  those  of  Ad  and 
At ,  respectively.  We  state  the  result  in  terms  of  the  /oo-norm  only  because  the  bounds  then 
have  a  simple  form. 

Remark  4.  We  have  computed  the  condition  numbers  k(V )  for  At  and  found  them  to  be 
large,  on  the  order  of  103  or  higher,  for  the  three  values  of  //,  with  streamline  upwinding  and 
h  =  1/16.  However,  the  presence  of  k(V )  in  these  bounds  is  an  artifact  of  the  proof  of  the 
Bauer-Fike  theorem;  there  are  more  subtle  analyses  ([8],  pp.  344ff),  as  well  as  bounds  that 
do  not  require  diagonalizable  matrices  [10].  We  have  observed  that  the  eigenvalues  of  At  are 
insensitive  to  perturbations,  and  we  believe  that  the  presence  of  k(V )  is  pessimistic.  This 
supposition  is  supported  by  the  experimental  results  described  below. 

To  demonstrate  that  inexact  preconditioning  is  effective,  we  consider  two  iterative  methods 
based  on  line-oriented  splittings  of  F.  The  first  uses  a  horizontal  line  Gauss-Seidel  relaxation: 
Let  F  =  H  —  R  denote  a  horizontal  line  Gauss-Seidel  splitting  of  the  block  convection- diffusion 
operator  F  derived  from  the  1-line  natural  left -to- right,  bottom-to-top  ordering  of  the  velocity 
grid.  Thus,  H  is  a  block  lower  triangular  matrix  consisting  of  the  block  diagonal  of  F  (a 
tridiagonal  matrix)  together  with  the  strict  block  lower  triangular  part  of  F.  (See  [17,  21]  for 
further  details.)  The  horizontal  line  Gauss-Seidel  method  for  Fw  =  v  performs  the  iteration 

wo  =  0,  wi+1  =  Wi  +  R-1(w  -  Fwi). 

For  k  steps  of  this  iteration,  the  approximating  matrix  is  F  =  F(I  —  (R_1i?)fc)_1 . 

It  has  been  observed  that  the  performance  of  relaxation  methods  of  this  type  can  be 
improved  if  the  sweep  direction  follows  the  underlying  direction  of  flow  [4].  Our  benchmark 
problem  has  a  circular  flow,  so  that  no  simple  line  relaxation  can  mimic  the  flow  direction 
throughout  0.  A  slightly  more  sophisticated  idea  is  to  use  an  alternating  line  relaxation.  For 
this,  let  F  =  V  —  T  denote  a  vertical  line  Gauss-Seidel  splitting  of  _F;  that  is,  if  P  is  a  permutation 
matrix  associated  with  the  mapping  from  the  natural  horizontal  line  ordering  of  grid  points  to 
the  natural  vertical  line  ordering,  then  PTVP  is  the  block  lower  triangular  part  of  PT FP.  One 
iteration  of  alternating  line  relaxation  consists  of  two  line  Gauss-Seidel  steps,  one  using  the 
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horizontal  splitting,  followed  by  one  using  the  vertical  splitting: 


wo  =  0,  wi+1/2  =  Wi  +  H  1(v-Fwi),  wi+1  =  wi+1/2  +  V  1(v  -  Fwi+1/2). 


i/=i 


v=l/10 


i/=l/100 


Fig.  6.  Performance  of  block  tridiagonal  inexact  preconditioners,  for  32  X  32  grid. 

Figure  6  shows  the  results  of  using  the  inexact  block  tridiagonal  preconditioners  with  both 
CfMR.ES(lO)  and  QMR,  to  solve  the  benchmark  problem  discretized  by  streamline  upwinding 
on  a  32  X  32  grid.  Results  for  inexact  block  diagonal  preconditioners  were  similar,  except  that, 
as  with  the  exact  preconditioners,  convergence  was  slower.  We  used  four  steps  of  horizontal 
line  relaxation  or  two  steps  of  alternating  line  relaxation,  so  that  both  inexact  preconditioners 
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perform  four  sweeps.  The  figure  also  shows  the  performance  of  the  exact  preconditioner,  whose 
cost  per  step  is  significantly  more  expensive.  For  example,  with  an  n  X  n  velocity  grid,  direct 
solution  using  a  bandsolver  requires  0(n 4)  operations,  whereas  each  inner  iteration  is  an  0(n2) 
computation.  We  see  that  the  use  of  inexact  preconditioners  in  place  of  the  exact  versions 
leads  to  little  degradation  of  performance  of  the  Krylov  subspace  methods.  For  example,  in  the 
convection-dominated  case  v  =  1/100,  QMR  with  alternating  line  relaxation  requires  roughly 
25%  more  iterations  than  with  the  exact  preconditioner.  For  the  diffusion  dominated  case  v  =  1, 
roughly  three  times  as  many  outer  iterations  are  required  with  the  inexact  preconditioners,  still 
leading  to  a  less  costly  computation.  Not  surprisingly,  alternating  relaxation  is  more  effective 
than  horizontal  relaxation,  especially  for  convection-dominated  problems.  We  remark  that  our 
goal  here  is  only  to  demonstrate  “proof-of-concept;”  many  other  techniques  for  approximating 
the  action  of  F~x  are  possible,  for  both  diffusion-dominated  and  convection-dominated  flow. 
See,  for  example,  [1,  4,  13,  15,  19]. 

In  conclusion,  in  this  paper  we  have  shown  that  appropriately  preconditioned  Krylov  sub¬ 
space  methods  can  be  used  to  solve  the  discrete  Oseen  equations,  a  linearized  version  of  the 
Navier- Stokes  equations,  with  convergence  rates  that  are  independent  of  the  mesh  size  used  in 
the  discretization.  This  work  generalizes  results  for  the  (self-adjoint)  Stokes  equations.  Precon¬ 
ditioning  techniques  derived  using  approximate  solution  methods  for  a  subproblem,  the  discrete 
convection- diffusion  equation,  have  modest  costs  and  lead  to  efficient  solution  algorithms  for 
the  linearized  Navier- Stokes  equations. 
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